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ABSTRACT 

We use a suite of cosmological hydrodynamic simulations including a self-consistent 
treatment for inhomogeneous reionisation to study the impact of galactic outflows and 
photoionisation heating on the volume-averaged recombination rate of the intergalactic 
medium (IGM). By incorporating an evolving ionising escape fraction and a treatment 
for self-shielding within Lyman limit systems, we have run the first simulations of 
"photon-starved" reionisation scenarios that simultaneously reproduce observations of 
the abundance of galaxies, the optical depth to electron scattering of cosmic microwave 
background photons r cs , and the effective optical depth to Lycc absorption at z = 5. 
We confirm that an ionising background reduces the clumping factor C by more than 
50% by smoothing moderately-overdense (A =1-100) regions. Meanwhile, outflows 
increase clumping only modestly. The clumping factor of ionised gas is much lower than 
the overall baryonic clumping factor because the most overdense gas is self-shielded. 
Photoionisation heating further suppresses recombinations if reionisation heats gas 
above the canonical 10,000 K. Accounting for both effects within our most realistic 
simulation, C rises from < 1 at z > 10 to 3.3 at z = 6. We show that incorporating 
temperature- and ionisation-corrected clumping factors into an analytical reionisation 
model reproduces the numerical simulation's r cs to within 10%. Finally, we explore 
how many ionising photons are absorbed during the process of heating filaments by 
considering the overall photon cost of reionisation in analytical models that assume 
that the IGM is heated at different rcdshifts. For reionisation redshifts of 9-10, cold 
filaments boost the reionisation photon budget by ~ 1 photon per hydrogen atom. 
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1 INTRODUCTION 

Much of the current interest in cosmological hydrogen reion- 
isation may be distilled to two related questions: (1) When 
did it occur?; and (2) What sources dominated the ion- 
ising photon budget? Theoretical efforts to address the 
first question have traditionally considered two observa- 
tional constraints: the optical depth to Thomson scatter- 
ing of cosmic microwave background (CM B) photons, r es , 
and the opacity of the Lyman-a forest |Trac fc Gnedinl 
120091 . and references therein). While these works have shed 
much light on the nature of inhomogeneous reionisation, 
their predictions have until recently been uncertain owing to 
the u nknown relation betw een ionising luminosity and halo 
mass (McQ umn et al.ll2007l ). With the recent installation of 
the Wide-Field Camera 3 on board the Hubble Space Tele- 
scope, a new generation of observations is now constrain- 
ing the abun dance and colors of ga l axies back to z = 8 
and beyond (|Finkelstein et all l20ld : iBouwens et al.l l201ll : 



Dunlop et al. 20121: Gonzalez et al. 201ll: Grazian et al.l 
201 it iMcLure et all l201ll : lOesch et all |2012| ). These ob- 
servations may be used to constrain the relationship be- 
tween star forma tion rate and halo mass via abundance - 
matching studies (|Trenti et alj|20ict iMunoz fc Loebll201ll ). 
Furthermore, the relationship between ionising luminosity 
and star formation rate is now constrained by the direct de- 
tection o f escaping Lyman continuum flux from high-redshift 
galax ies (|Shaplev et aHl200d ; ls"iana et alj2010tlNestor et all 
120111 ). Hence there are now tenuous observational links con- 
necting halo mass with ionising luminosity, a major step to- 
ward understanding the relationship between high-redshift 
galaxies and their environment. 

By combining observations of galaxies, the CMB, 
and the Lyman-a forest, it is possible to ask whether 
the observed ionising sources and reionisation h istory 
are compatible (for example , lHaardt fc Madaul 120121 : 



iKuhlen fc Faucher-Giguerell2012t ). Two uncertainties hinder 
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Figure 1. Our ionisation- and temperature-corrected IGM 
clumping factor at z = 6 in comparison to previous calculations. 
Our work indicates that the clumping factor lies within the range 
2.7—3.3, which is consistent with the low values that have been 
found in recent works. 



progress toward answering this question. The first remains 
the fraction of ionising photons / osc that escape into the 
IGM: direct constraints are only available at lower redshifts 
and higher luminosities than are relevant to reionisation. 
Furthermore, / csc is sensitive to parsec-scale physical pro- 
cesses that occur deep within galax ies' interstellar media and 
is difficult to model numerically (jFernandez fc Shulll 1201 J . 
and references therein). The second uncertainty involves the 
recombination rate of the reionisation-epoch IGM. The un- 
known ratio of the IGM's true recombination rate to its 
hypothetical rate under the assumption of uniform density 
and temperature is often referred to as the IGM clumping 
factor C. We will define C explicitly in Section [3] for the 
present the important point is that it is proportional to the 
volume-averaged hydrogen recombination rate in the IGM. 
High values of C require abundant faint galaxies and/or a 
large / csc while low values of C are easier to reconcile with 
existing observations. C is sensitive to the IGM's density, 
ionisation, and temperature fields, hence an accurate esti- 
mate requires three-dimensional numerical simulations. 

To this end, C has been studied numerous times us- 
ing numerical simulations. We compile predictions of the 
clumping factor at z — 6 from a variety of works in Fig- 
ure Q] The early work of iGnedin fc OstrikeJ (|l997h yielded 
a clumping factor of 30-50 (their Figure 2) depending on the 
choice of definition. However, this clumping factor included 
the recombinations that occur in gas whose ionisation state 
is dominated by the local field rather than by the global ex- 
tragalactic ultraviolet ionising background (EUVB), hence 
it is an overestimate of the recombination rate within the 
gas that is conve ntio nally associated with the IGM (see 
also lGnedinll2OO0l and|M cQuinn et al.|[201ll ). For the same 



reason, iTrac fc Cenl (|2007h also inferred a clumping factor of 
roughly 30 from simulations that subtended a significantly 
larg er cosmological v olume . 

Illiev et~al] (|2007l ) applied this insight to calculate an im- 
proved clumping factor using high-resolution N-body simu- 
lations whose spatial resolution was comparable to the Jeans 
length in the IGM. By excluding matter within virialized ha- 
los, they derived a redshift-dependent fitting function that 
climbs to ~ 10 at z = 6, a val ue which was also found i n 
the recent N-body simulations of lRaicevic" fc Theunsi(|201lf ). 
This work confirmed that physically-motivated density cuts 
reduce C. However, as the authors noted, their calculation 
did not account for the fact that photoionisation heating 
(photoheating) tends to smooth the diffuse IGM. We will 
confirm in Figure [5] that this leads to an overestimated 
clu mping factor. 

IPawlik et ail (|2009r ) corrected this problem by running 
cosmological hydrodynamic simulations that model reionisa- 
tion using a spatially-homogeneous EUVB. They used sev- 
eral plausible density cuts to calculate the IGM clumping 
factor as a function of redshift. We show their predictions 
assuming that the IGM corresponds to overdensity cuts of 
A = p/{p) < 50 (lower range) and A < 100 (upper range). 
This work generally confirmed the impression that consider- 
ing only low-density gas and accounting for the tendency for 
the photoheating of diffuse regions both reduce C. However, 
it suffered from two drawbacks. First, their use of a spatially- 
homogeneous EUVB pressurized regions that should in real- 
ity be self-shielded, leading to uncertainty in the gas density 
distribution from which the predicted clumping factors were 
derived. Second, their analysis did not account explicitly for 
the IGM temperature, which could suppress the recombina- 
tion rate further. 

In a qualitatively similar work, IShull et al.l (|2012l ) ran 
a cosmological hydrodynamic simulation with a spatially- 
homogeneous EUVB and computed the temperature- 
corrected clumping factor within gas that satisfied cuts on 
baryon overdensity, temperature, metallicity, and ionisation 
state . Although th e se ext ra cuts render direct comparisons 
with IPawlik et all (|2009l ) uncertain, the inferred range is 
consi stent. Both works are also in reasonable agreement 
with iFaucher-Giguere et al.l (l200ct). who used the gas den- 
sity distribution of Miralda-Escude et al.l (|2000h to find that 
C =3-4 if the clumping factor is averaged only over gas with 
overdensity A less than 67. 

The uncertainty in the derived clumping factor may be 
reduced by using radiation transport simulations to identify 
directly the gas that is self-shielded. iMcQ umn et all (|2011l ) 
use this approach in post-processing to derive a clumping 
factor of 2.4-2.9. Their calculations do not account f or the 
IG M temperatu r e and, like the works of lPawlik et all (|2009h 
and IShull et al.l (|2012T l. begin with hydrodynamic simula- 
tions in which the moderately-overdense gas is incorrectly 
exposed to the EUVB (although they confirm directly that 
shielding gas that is dense enough to form stars has no im- 
pact). Nonetheless, their results are clearly consistent with 
a relatively low clumping factor. 

The goal of the present work is to combine all of these 
ideas into a single radiation hydrodynamic framework that 
accounts realistically for the time and spatial dependence of 
the IGM's density, ionisation, and temperature fields. Our 
fiducial simulation reproduces the observed UV luminosity 
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function of galaxies at z ^ 6 (|Finlator et a.1.1 [201 if ) and is 
marginally consistent with the effective optical depth to ab- 
sorption of Lya photons and r os . Hence, our estimated IGM 
recombination rate is constrained by a variety of observa- 
tions. As Figure [1] indicates, we will confirm that effects 
that have been neglected in recent works based on hydro- 
dynamic simulations such as incorrectly heating dense re- 
gions and ignoring the gas temperature are indeed small. Of 
course, the purpose of clumping factors is to inform analyt i- 
cal models of reionisation (for example. iMadau et al.|[T999t ). 
We will therefore verify that analytical reionisation calcula- 
tions incorporating clumping factors do, in fact, reproduce 
the behavior of numerical simulations. This check will as- 
sess directly the magnitude of the systematic uncertainties 
inherent in calculations that ignore higher-order effects such 
as shadowing, source clustering, and the light-travel time of 
ionising photons. 

Having verified that our analytical model reproduces 
the behavior of our numerical model reasonably well, we will 
use it to ask how many photons are absorbed by filamentary 
overdensities in the diffuse IGM. This will quantify the error 
introduced in models that compute reionisation using IGM 
density distributions in which the gas is effectively preheated 
prior to the onset of reionisation. 

The outline of our paper is as follows: We begin by in- 
troducing our suite of simulations in Section [21 We compare 
different definitions of the clumping factor in Section [3] and 
explore how feedback effects modulate their values in Sec- 
tion U In Section [SJ we incorporate our derived clumping 
factors into an analytical model for reionisation and test how 
well analytical models of reionisation reproduce the behav- 
ior of numerical calculations. We apply this model to study 
how many photons are absorbed in filaments in Section [6] 
Finally, we discuss our results in Section [7] and summarize 
in Section [5] We introduce two improvements to our treat- 
ment of the Eddington tensors in Appendix [X] and quantify 
resolution limitations in Appendix FBI 



2 SIMULATIONS 

2.1 Star Formation and Outflows 

Our numeric al methods, input p hysics, and cosmology 
are similar to iFinlator et al.l (|201 lh with the exception of 
our adopted treatments for an evolving / esc , subgrid self- 
shielding (see below), and two numerical optimizations re- 
garding the treatment of the Eddington tensor field (Ap- 
pendix fS}. We model the growth of structure and the feed- 
back processes that couple galaxies with the IGM using our 
custom vers ion of the cosm ological galaxy formation code 
Gadget-2 (|Springell I2005T ) . Gadget-2 implements a for- 
mulation of smoothed particle hydrodynamics (SPH) that 
simultaneously conserves entropy and energy and solves 
for the gravitational potential with a tree-particle-mesh al- 
gorithm. Dense gas cools radia tively using the p rimordial 
cooling processes in Table 1 oflKatz et all (| 19961) a nd th e 
metal-line cooling tables of ISutherland fc Dopital {l993h . 
whic h assume co l lisiona l ionisation equilibrium. In contrast 
with iKatz et al.l (|l996p . however, we evolve the ionisation 
states of hydrogen and helium and the electron abundance 
simultaneously with the cooling equations using a nested 



subcycling approach w hose timestep is limite d by a chemical 
Courant condition; see lFinlator et all (|201lf ) for details. We 
initialize the IGM temperature and neutral hydrogen frac- 
tion to the values appropriate for each simulation's initial 
redshift as computed by recfast (|Wong et al.l 120081 ) . and 
we assume that helium is initially compl etely neutral. We 
genera te the initial density field using an lEisenstein fc Hul 
jl999l ) power spectrum at redshifts of 249 and 319 for simu- 
lations subtending 6 and 3 /i _1 Mpc, respectively. All simu- 
lations assume a cosmology in which Q.m = 0.28, Qa = 0.72, 
Q b = 0.046, h = 0.7, a 8 = 0.82, and the index of the pri- 
mordial power spectrum n = 0.96. 

Our goal of modeling galaxies and reionisation simul- 
taneously requires a treatment for the ability of cool gas 
to form stars. We ado pt the subgrid two-ph a se int erstellar 
medium treatment of ISpringel fc Hernquistl (|2003l ). which 
can be tuned to reproduce the observed relatio n between 
the s urface densities of gas and star formation (|Kennicuttl 
Il998l ). The physical density threshold for star formation is 
0.13 cm -3 . This value is motivated by observations of a crit- 
ical density for the onset of star formation and lies within 
the range at which the th ermo- gravitat ional instability is ex- 
pected to become active (Schavc 2004). Varying it has only 
a mi nor impact on the predicted star formation rate den- 
sity {Schave et alJlioiol ). We account for metal enrichment 
owing to supernovae of Type s II and la as well as a symp - 
totic giant branch stars; see lOppenheimer fc Pavel {2008) 
for details. 

Allowing cold gas to form stars without any feed- 
back invariably results in overproducing th e observed 
reion i sation-epoch star fo rmation rate density {Dave et al.l 
120061 ; IFinlator et all [20111 ). The accepted solution is to al- 
low feedback from massive stars to expel star-forming gas 
from galaxies. As the spatial resolution necessary to form 
such outflows self-consistently i s beyond the reach of cur- 
rent cosmolog i cal s i mulations dMac Low fc Ferraral 1 19991 ; 
iHopkins et ail |2012| ; iPowell et all l201ll ). we impose star 
formation-driven outflows by stochastically applying kicks 
to star-forming gas particles. The amount of gas kicked 
per unit stellar mass formed and the kick velocities are ad- 
justed to reprod uce the expected sca lings from momentum- 
driven outflows {Murray et al.l [20051 ). We temporarily dis- 
able hydrodynamic forces in outflowing gas in order to mimic 
the way in which outfl ows escape through holes i n higher- 
resolution simulations {Mac Low fc Ferraral 119991 1. Hydro- 
dynamic forces are restored once the gas has traveled for 
1.95 x 10 4 km s -1 /u million years (where v is the kick veloc- 
ity) or its density drops below 10% of the threshold density 
for star formation. In practice, this allows outflowing gas to 
reach galactocentric radii of 50-100 physical kpc and then 
re-ac crete onto the central galax y on a timescale of ~ 1 
Gyr {Oppenheimer fc Pavel l20ul ). The distance that out- 
flows travel before re-accreting is large compared to the typ- 
ical radius at which hydrodynamic forces are enabled, hence 
the decoupling prescription has little impact on outflows 
once they escape the ISM. The long re-accretion timescale 
enables expelled gas to modify the IGM recombination rate 
and opacity depending on the density of outflowing gas; ex- 
ploring these possibilities is one of the goals of the present 
work. 
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2.2 Radiation Transport 

We evolve the EUVB on a fixed Cartesian grid using the 
moments of the cosmological radiative transfer equation. We 
close the moment hierarchy with a variable Eddington ten- 
sor field that is updated periodically via a time-independent 
ray-casting calculation. This approach accounts accurately 
for the inhomogeneous o pacity field as long as the updates 
occur frequently enough (jFinlator et al.ll2009al ). Each cell's 
ionising luminosity is given by a sum over the luminosities 
of its star-forming gas particles. Their luminosities, in turn, 
depend on their star formation rates and metallicities, where 
the metallicity depe ndence follows an analytical fit to Table 
4 of lSchaererl (|2003l ). 

The fraction of ionising photons that escapes into the 
IGM in order to participate in reionisation / osc is unknown. 
Recent work suggests that it must vary steeply with red- 
shift in order to bring observations of galaxies, the Lyman- 
q forest, and the C MB optical depth to Thomson scatter- 
ing into agreement (llnoue et al.l 20061; Finlator et alJ 201 it 
lHaardt fc Madaul |2012| ; IKuhlen fc Faucher-Giguerel (20121 ). 
We adopt a strongly redshift-dependent f e . 
lations: 



in our simu- 



/csc — 



/gsc,5 ( 1 g Z ) 

1.0 



2 < 10 
2 > 10 



(1) 



Here, the normalization / csc ,5 sets the escape fraction at 
z = 5 and k controls how strongly / csc varies with red- 
shift. For simulations with outflows, we use / csc ,5 = 0.0519 
and k — 4.8. The normalization i s consistent with ob- 
servations of Lyman break galaxies (| Nest or et al.|[20lH ) at 
2 ~ 3, while the power-law index lies w ithin the range (1-6) 
that IKuhlen fc Faucher-Giguerel (|2012f ) indicate is required 
by observations if the UV luminosity function of galaxies 
extends to absolute magnitudes of -15 (as required by ob- 
servations of gamma-ray bursts; llrenti et af] 12012; ) . With 
this parametrization, the effective optical depth to Lyman- 
a absorption at z = 5 is r a =3.1, marginal ly consistent with 
the observed range of 2-3 (|Fan et alj|200l v Meanwhile, the 
predicted ionising emissivities are (4.6 , 7.0) x 10 50 s -1 Mpc -3 
at 2 =(5, 6). At z = 5, this is consistent with the ob - 
served range of 4.3 ± 2.6 (|Kuhlen fc Faucher-Giguirell2012l ). 
At z — 6, however, it overshoots the observed limit < 2.6, 
suggesting that the predicted emissivity should strengthen 
rather than declining toward the end of reionisation. Note 
that our model is not unique in failing to reproduce the 
observed emissivities at 2 — 5 and z — 6 simultane- 
ously; this evolution is evidently quite strong compared 
to expectations from models ( see, fo r example, Figure 10 
of IKuhlen fc Faucher-Giguerel I20l3 or Figures 8 and 15 
of lHaardt fc Madau 2012). The resulting reionisation his- 
tory corresponds to an optical depth to Thomson scattering 
of 0.071, only slightly below t he observed la confid ence in- 
tervals of r cs = 0.088 ± 0.015 jKomatsu et alj|201ll ). 

Simulations without outflows require a steeper power- 
law index and a lower normalization because they ove rpro- 
duce the observed galaxy abundance |Dave et al. lixi). We 
adopt /csc, 5 = 0.0126 and k = 7.21. With this dependence, 
they predict emissivities of 5.4 and 7.8 x 10 50 s _1 Mpc~ 3 at 
z = 5 and 6, respectively, once again indicating good agree- 
ment at z = 5 but incorrect evolution to higher redshift. 
The predicted Lyman-a optical depth at z = 5 is 2.2, well 



within the observed range. The predicted Thomson scatter- 
ing optical depth is 0.074, within the observed lcr confidence 
intervals. 

The physical interpretation of an evolving / csc could 
take many forms. For example, direct measurements of ion- 
ising continuum flux from galaxies at z as 3 suggest that cur- 
rent stellar population models underestimate t he ionising lu- 
mino sity of low-metallicity Population II stars (|Nestor et al.l 
l201lT ). If confirmed, this observation could imply an evo- 
lutionary trend to higher ionising luminosity at high red- 
shift that mimics an evolving / csc . Alternatively, observa- 
tions indicate that gala xies grow bluer and by i nference less 
dusty beyond 2 = 4 |Finkelstein et all |2010| ) ; this could 
drive u p f csc if dust domi nates the absorbing column (al- 
though iGnedin et al. I l2008l argue that / csc is dominated by 
the ISM's geometry rather than its dust content). Finally, 
galactic outflows can fill a halo with gas, potentially creat- 
ing a significant optical depth to ionising photons. If these 
screens grow weaker with increasing lookback time, then the 
fraction of ionising photons that travel unimpeded from the 
edge of a galaxy's ISM to the virial radius could grow accord- 
ingly; this would mimic an evolving / esc . We defer detailed 
investigation of these possibilities to future work; for the 
present, we simply choose a parametrization that will im- 
prove the realism of our simulated density field. We discuss 
drawbacks to this approach in Section [7] 

We assume that each photoionisation deposits 4.08 eV 
of latent heat into the IGM. This heats newly-reionised gas 
to ~ 15, 000 K, which lies within the range that is expected 
for a Population II stellar spectrum. Each radiation trans- 
port cell's opacity is given by its volume-averaged neutral 
hydrogen density (counting only the gas that is not self- 
shielded; see below) multiplied by the cross-section for ioni- 
sation at an energy of 17.68 eV. 

2.3 Self-shielding Within Overdense Systems 

Two central challenges in using numerical models to com- 
pute the IGM recombination rate involve deciding which gas 
belongs to the IGM and ensuring that overdense gas is al- 
lowed to self-shield against the EUVB. 

The need to define the IGM in a physically-motivated 
way is illustrated by the results of the pathbreaking radia- 
tion hydrodynamic simulations of lGnedin fc Ostriker] (|l997h 
and iGnedinl (|2000h . These simulations faithfully modelled 
the growth of ionised regions as well as the hydrodynamic 
response of the heated gas. However, the volume-averaged 
recombination rates were computed over all the gas in the 
simulation, leading to recombination rates that were much 
larger than what is plausible for the moderately overdense 
(p/(p) < 100) gas that is conventionally associated with the 
IGM. It also led to a strong dependence of the recombination 
rate (and the reionisation photon budget) on the simulation 
dynamic range, w ith higher-res olution simulations absorb- 
ing more photons (|Gnedinll2000l ). 

Improvement results from using simple density cuts to 
isolate low-density gas in post-processing from precomputed 
simulations. This approach has been applied both to N- body 
simulations (|lliev et al.ll2007l ; iRaicevic fc Theunsl 2011) and 
to cosmological hydro dynamic simulations ( Pawlik et al.l 
120091 ; IShull et alj[20l3 ). It reduces C as expected, but the 
results remain dependent on the precise cuts used to de- 
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fine the IGM. In the case of hydrodynamic simulations, 
additional uncertainty results from assuming a homoge- 
neous EUVB, which may pressurize regions that are in re- 
ality dense enough to self-shield and does not account for 
the spatial inhomogeneity of reionisation. Both of these 
effects could modify the IG M's density and temperature 
fields (jFurlanetto fc Ohll2009h . which in turn impact the re- 
combination rate. 

In this work, we combine these two approaches by us- 
ing radiation hydrodynamic simulations of the reionisation 
epoch that incorporate a physically-motivated treatment for 
self-shielding within overdense regions. Effectively, we in- 
corporate directly into our simulations the definition that 
the IGM consists only of that gas whose ionisation state is 
not determined by local sources; denser systems are allowed 
to remain neutral. By isolating overdense systems from the 
EUVB, we (1) improve the accuracy of our predicted den- 
sity and temperature fields; (2) prevent our simulations from 
overestimating the reionisation photon budget; and (3) obvi- 
ate the need to choose a self-shielding threshold overdensity 
(such as 50 or 100) when deriving the IGM recombination 
rate in post-processing: the threshold evolves on-the-fly in 
a way that follows the (local) gas temperature and (global) 
EUVB amplitude. 

We model the self-shielding of overdense regions us- 
ing a subgrid recipe because our radiation transport solver 
does not resolve them spatially even though our hy- 
drodynamic solver does. For example, Lyman limit sys- 
tems (with neutral column densities of ~ 1 17 cm~ 2 ) have 
a characteristic size of 10 physical kpc (|Schavd l200ll ; 
iMcQ uinn et al.ll201 lh and must be treated with 5-10 resolu- 
tion elements ( Gnedin fc Fanl2006l ') . In our higher-resolution 
simulation (r6n256wWwRT32), the gravitational softening 
length is 0.1 kpc (Plummer equivalent) while the radia- 
tion transport solver achieves a resolution of 38 kpc (where 
both numbers are in physical units at z = 6). Hence while 
our hydrodynamic solver certainly resolves Lyman limit sys- 
tems, our radiation solver's resolution is a factor of w 20 too 
coarse. 

We overcom e this limitation by following the argument 
of ISchavel (|200ll ). This work showed that, if Lyman limit 
systems are in hydrostatic equilibrium, then they can be 
iden tified directly with ga s above a critical overdensity (see 
also iMcQuinn et al.ll2011T ). The threshold, An s , is given by: 



A lls = 25 



f— Y 

V10 4 K J 



l + z 



2 x 10" 



C2) 



at the Lyman limit. Here we have assumed that hydrogen 
is fully ionised and heli um is singly-ionis e d, as expected 
for so ft ionising sources (|Ciardi et all 120121 ; iFriedrich et al.l 
120121 '). Fhi represents the volume- averaged EUVB and does 
not include the influence of local sources. This estimate is 
within a factor of two of what others have found (for ex- 
ample, iBolton fc Haehneltll2007l '). We do not expect slight 
differences owing, for example, to our assumed recombina- 
tion rate to affect our results. We enforce An s 2 in order to 
preserve self-shielding at z > 10; this floor is not expected 
to affect our results. We use this approximation to model 
subgrid self-shielding as follows: For each gas particle whose 
density lies below the physical threshold for star formation 
(0.13 cm -3 ), we compute An s using the local gas tempera- 
ture, the case-A recombination rates, and the global mean 



ionisation rate per hydrogen Fhi. We then assume that the 
optical depth to the diffuse background varies as a power-law 
of the overdensity 



Tr = 





To 



An, 



A < An 
A > A„ 



(3) 



In the optically-thin limit, To = 1 and b = 1.5 jSchaydbOOll ). 
In reality, rr increases more rapidly with density because 
the neutral fraction in the partially shielded outskirts of an 
absorber is larger than in the optically-thin case. We have 
experimented with both b = 1.5 and b = 3 and found little 
difference in practice. We adopt the fiducial choices to = 1.0 
and b = 3 in this work. 

We implement Equations [2}[3] directly into our simu- 
lations so that overdense gas is shielded from the EUVB 
on-the-fly. Our use of different discretizations for the radia- 
tion and hydrodynamic solvers requires that we do so in two 
ways. When updating the ionisation state of a given parti- 
cle, we attenuate the radiation field of its host cell by the 
factor e _Tr , where rr is evaluated using the global EUVB 
and the particle's temperature. Conversely, when assembling 
the opacity field on the radiation solver's grid, we reduce 
each particle's contribution to its host cell's opacity by the 
same factor. We include the opacity of each gas particle's 
self-shielded region by treating it as an opaque sphere with 
volume equal to /shi(1 — e _Tr )AU, where AV is the ra- 
tio of the gas particle's mass to its density, xni is its neu- 
tral fraction, and / is a parameter that is tuned via high- 
resolution calculations to 1/8. For example, if e~ Tr = 0.5 for 
a particle, then its ionisation rate is 50% of that of its host 
cell's, and only 50% of its neutral hydrogens contribute to 
the cell's opacity. We also augment the cell's opacity by the 
amount 4.84(0.5:eh iAV)~ 2/3 /A cclh where A cen is the cell's 
area. Note that the opacity from partially-neutral gas parti- 
cles whose overdensity falls below the self-shielding thresh- 
old is included directly in the radiation transport solver as 
in our previous work. 

Intuitively, this treatment divides self-shielded regions 
into an optically thin skin and an optically-thick core. 
The overdensity range over which a region transitions from 
optically-thin to optically-thick follows from Equation [3] 
The opacity owing to the optically-thin region is distributed 
uniformly throughout the radiation transport cell while the 
core is treated as a photon sink. This approximation pre- 
serves the ability of dense gas to remain self-shielded within 
a coarse radiative transfer grid. It will enable radiative hy- 
drodynamic simulations to subtend cosmological volumes 
even as they treat the radiation field's small scale struc- 
ture with sufficient detail to model observables such as the 
abundance of low-ionisation metal absorbers and the post- 
reionisation IGM opacity. 

In reality, regions that are dense enough for star forma- 
tion to occur (A > 1000) contain a mixture of ionised and 
neutral gas owing to the local radiation field from massive 
stars. They remain entirely neutral in our simulations ow- 
ing to our self-shielding prescription. Accounting for their 
reco mbinations would lead to much larger clumping fac- 
tors (|Gnedin fc Ostrikerlll997l : lGnedinll200oh . However, the 
corresponding absorptions should be identified with an ion- 
ising escape fraction / csc that is less than unity rather than 
an increased clumping factor because they are not domi- 
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Figure 2. The threshold overdensity for self-shielding An s as a 
function of redshift within our fiducial simulation at a charac- 
teristic temperature of 10 4 K. An s peaks at over 200 during the 
heyday of rcionisation before shrinking to roughly 150 by z = 5 
owing to declining escape fractions. 



nated by the EUVB. We encapsulate their impact via an 
assumed / csc ^ 1 that is tuned using observational insight 
(see above) and defer the challenging task of modeling / esc 
directly to future work. 

Figure [5] shows how An s evolves in our fiducial simula- 
tion for a characteristic temperature of 10 4 K. Note that this 
figure is illustrative; in reality our simulations compute this 
threshold on-the-fly using the local temperature. An s grows 
from values near unity at the onset of reionisation to values 
near 200 by the time that reionisation has completed. Gas 
that is denser than this threshold sees an attenuated EUVB, 
suppressing both the reionisation photon budget and the 
volume-averaged recombination rate. 

In detail, the values in Figure [2] may be somew hat over- 
estimated (see, for example. iMcQuinn et al.ll201ll ). This oc- 
curs because our monochromatic radiation transport solver 
uses a smaller photoionisation cross section a than the value 
at the Lyman limit ctll in order to match the assumed pho- 
toheating rate of 4.08 eV per photoionisation. In particu- 
lar, Aii s vari es with the p hotoionisation cross section a as 
(cr/Wr 273 llSchaydbOOfj ). hence our simulations may over- 
estimate Aii s by a factor of « 1.7. This is ro ughly the differ- 
ence b etween our value at z = 6 and that of lMcQuinn et all 
|201ll ). The overall impact of this limitation on our inferred 
value of C is weak; we will quantify it in Section [7] 

Table [T] introduces our simulation suite. The simulation 
names encode the physics and numerical resolutions. For 
example, the r6n256wWwRT16 simulates a 6/i _1 Mpc vol- 
ume (r6) with 2 x 256 3 particles (n256), includes outflows 
(wW) and evolves the EUVB on a Cartesian grid with 16 3 
cells (wRT16). The first four simulations explore the sen- 
sitivity of gas clumping to photoheating and outflows. The 



next three explore convergence with respect to our radiation 
and gas solvers as well as the cosmological volume. The final 
simulation assumes an optically-thin EUVB for comparison 
with previous work. Throughout this work, we will use the 
r6n256wWwRT16 simulation as our fiducial case. 

Our simulations allow us to explore how the IGM's tem- 
perature impacts its recombination rate within the context 
of a model that treats photoheating from Population II stars 
self-consistently (Section 14. ip . As a check on whether our 
predictions are realistic, we compare the temperature at 
mean density to recent observations. iBolton et alj (|2012l 'l 
used measurements of the Doppler widths of Lycv absorption 
lines along the line of sight to seven quasars to infer that the 
IGM temperature at mean density is log(T) = 3.85 ± 0.08 
at z = 6.08 ± 0.33 (their Table 3, fiducial model). Our fidu- 
cial simulation predicts log(T, z = 6) = 3.93 ± 0.10. Here, 
we report the median log(T) over particles with overden- 
sity —0.2 < log(A) < 0.2; switching from the median to the 
mean changes results by ~ 0.1% because the scatter in the 
temperature distribution is not large. This comparison is of 
course incomplete because it does not account for a variety 
of observational systematics. Nonetheless, the fact that the 
simulation reproduces the observed value within the errors 
indicates that the predicted temperature fields are plausible. 



3 WHAT IS GAS CLUMPING? 

The motivation for defining clumping factors comes from a 
desire to take the volume-weighted mean of the ionisation 
rate equations. For example, the rate of change of the neutral 
hydrogen abundance nm is given by 



dt 



— rmiiHi + fc 2 nHiin — fcinmne, 



(4) 



where Phi is the photoionisation rate per hydrogen atom, 
ki is the collisional ionisation rate, and ki is the (case A) 
recombination rate. Each term on the right side is nonlinear, 
hence, in principle, it is not possible to compute their spa- 
tial averages unless the cross-correlation of the abundances 
of nm and nmi with each other as well as with the radia- 
tion and temperature fields are known (since k\ and &2 are 
temperature-dependent). It is convenient to encapsulate this 
subgrid information with clumping factors. For example, the 
spatially averaged recombination rate is approximated as 



(feriHimo) = CHII,T b (n.Hii)(ne)(fc 2 ), 



(5) 



where angle-brackets indicate averages over the entire simu- 
lation volume. Similar clumping factors could be defined for 
each of the terms on the right-hand side of Equation U In 
practice, however, analytical calculations bypass the need 
for the first term's clumping factor by assuming that the 
volume-averaged Phi equals the ionising emissivity divided 
by the mean hydrogen abundance. Meanwhile, the third 
term can be neglected because Thi 3> kin e except within 
shocks. This means that, along with the ionising escape frac- 
tion, Ch ii, T b is one of the major uncertainties that hamper 
efforts to connect the observed abundance of ionising sources 
(such as Lym an Break galaxies) with the ionisation state of 
the IGM (see lKohler et al1l2007l for an expanded discussion 
of the other clumping factors). In Section [5j we will show 
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Figure 3. The evolution of four different clumping factors in our 
fiducial simulation volume. The clumping factors of all gas whose 
density exceeds the star formation threshold reaches 15 by z = 6. 
Accounting for the fact that overdense gas is self-shielded and 
that ionised gas is somewhat warmer than 10 4 K leads to much 
lower clumping factors (see text for discussion). 

that the numerical results may be reproduced by an analyt- 
ical model that considers only the recombination clumping 
factor even though the full numerical calculation accounts 
for all three terms. 

Equation [S] shows that the clumping factor depends 
on the density, ionisation, and temperature fields. Partial 
clumping factors are often used when one or more of these 
fields is unavailable. In Figure [3] we show how the choice of 
definition impacts the clumping factor's value in our fiducial 
simulation. The black solid curve shows the evolution of the 
clumping factor averaged over all gas whose density lies be- 
low the adopted threshold for star formation (0.13 cm -3 ), 
Cb = (p 2 )/{p) 2 - It increases monotonically as overdense re- 
gions collapse, with a brief enhancement around the reion- 
isation redshift. The enhancement occurs when the EUVB 
grows stron g enough to photoeva porate gas from photosen- 
sitive halos ( Finlator et alj|201ll ). temporarily boosting the 
mass of gas whose density lies just below the threshold for 
star formation. The clumping factor of all gas (including 
star forming gas) decreases at this epoch as expected, but 
the clumping factor of gas that lies outside of galaxies — 
that is, Cb — increases until the newly-pressurized gas has 
expanded into hydrostatic equilibrium. 

The clumping factor within ionised regions — which is 
the important quantity for reionisation — is expected to be 
different from Cb because reionisati on begins in overdense 
regions and progress es into voids (|Furlanetto et all |2004| ; 
iMcQuinn et al.ll2007t ) in such a way that the last regions t o 
remain neutral are overdense (jMiralda-Escude et al.ll2000h . 
Furthermore, gas that is self-shielded does not contribute to 
recombinations. To explore these factors, we compute the 



clumping factor of electrons and protons 

Chii = (n c n H ii)/(n ){nHii). (6) 

Ch ii can be computed in either of two ways. One way 
is to compute the averages over all ionised gas. Unfortu- 
nately, this introduces uncertainty owing to the choice of 
ionisation threshold. For example, averaging only over gas 
particles that are more than 50% ionised gives a slightly 
different answer than averaging only over particles that are 
more than 95% ionised because the simulated IGM is not 
a strictly two-phase medium. The alternative is to take av- 
erages over the entire simulation volume and then weight 
by the ionised volume fraction xnu.v- This can be under- 
stood by recalling that, in a partially-ionised universe with a 
uniform density, Chii = 1/:ehii,v- Therefore, xhii,vChii in- 
dicates the recombination rate in units of the recombination 
rate of a homogeneous universe that has the same ionisation 
state, which is the goal of the clumping factor. 

We show this quantity with a red dotted line. The 
tendency for overdense regions to remain neutral sup- 
presses £hii,vChii below Cb at all times, with reionisation- 
epoch values that lie below ~ 6. This difference immedi- 
ately reveals the importance of self-shielding, which cre- 
ates a boundary between the ionised IGM and the neu- 
tral or locally-ionised regions near halos. It can be mod- 
eled directly by spatially resolving physical scales of ~ 1 
kpc (jGnedin fc Fan! 120061 ). However, this incurs significant 
computational expense, limiting the simulation's cosmolog- 
ical volume. Our approach of coupling a somewhat coarse 
grid for the radiation solver with a subgrid prescription for 
self-shielding may open up the possibility of simulating cos- 
mological reionisation within large volumes while treating 
the IGM's thermal and ionisation states faithfully. 

Chii is an improved description of the volume-averaged 
recombination rate over Cb, but it remains incomplete 
because the recombination rate also depends on the 
gas temperature. To illustrate this, we use a blue dashed 
curve to show the temperature-corrected clumping factor 
^h ii.vChii.Tk , comp uted following Equation (referred to 
bv lShull et al.|[2012l as Crr). This curve remains near zero 
until reionisation is well under way because the recombi- 
nation rate per proton is spatially anti-correlated with the 
abundance of ionised gas. Once 2hii,v exceeds « 0.5, it be- 
gins to climb because (fe) is no longer suppressed by cold 
neutral regions. However, it does not reach Chii owing to the 
lingering presence of self-shielded regions. Following reion- 
isation, a;Hii,vCHii,T b climbs slowly because the volume- 
weighted mean temperature is supported by the slow photo- 
heating of filaments, which suppresses (fe). It never exceeds 
3, reflecting the tendency of photoheating to suppress re- 
combination rates. 

£Hii,vCHii,T b provides a reasonable description of the 
clumping factor within ionised regions once reionisation is 
well under way, but it is less informative at early times 
(a;Hii,v < 0.1) when the high recombination rate in the pre- 
dominantly cold, neutral IGM suppresses CHii,T b despite 
the high clumping factor of ionised regions. We may com- 
pute a more informative clumping factor by re-evaluating 
CHii,T b using only those regions (that is, SPH particles) 
whose ionised mass fraction is greater than 0.95. This clump- 
ing factor (green dot-dashed) is somewhat larger at earlier 
times because it is not suppressed by the high recombation 
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rates per proton in regions that have not yet been heated. It 
merges with a5Hii,vCHii,T b once £hii,v ~ 1. A drawback of 
CHii,T b is that its value depends on the threshold ionisation 
fraction that is used to compute it. 

While CHii,T b emerges naturally from a spatial average 
of Equation 2] it is difficult to compare with observations 
because the mean IGM temperature (and hence (£2)) at 
2 — 6 remains poorly-constrained; this is of course equally 
true even if CHii,T b is computed only over regions whose 
ionisation fraction exceeds some threshold. To resolve these 
problems, we also compute the "observational temperature- 
corrected clumping factor of ionised gas" 

n (n e nHiifc 2 (r)) 

°HII,10 4 K = 1 \j n — TTTTTTTT \' ) 

(n e )(n H 11) fa (10 4 K) 

in which we replace the volume-averaged recombination rate 
(fa) in the denominator with the recom bination rate at 
10 4 K. This clumping factor (referred to bv lGnedTnlbOOOl as 
Chii) expresses the mean recombination rate without requir- 
ing knowledge of the topology of reionisation or the temper- 
ature of ionised gas. It is slightly higher than CHii,T b follow- 
ing reionisation, reflecting the tendency for gas to cool below 
10 4 K once reionisation is complete (that is, (fa) > fe(10 4 K) 
for 2 < 6). The function 

a&m,vC H n,io*K = 9-25 - 7.21 log 10 (l + z) (8) 

fits its evolution to within 30% for 2 ^ 15 and 10% for 
2 ^ 10. 

Note that uncertainties in the recombination rate trans- 
late directly into uncertainties in the clumping factor. For 
example, our adopte d (case A) recombinat ion rate is ~ 20% 
higher than that of iHui fc Gnedinl j 19971 ) in the tempera- 
ture range 10 3 -10 5 K. This means that ionisation fronts do 
not penetrate as far into overdense regions in our simula - 
tion as they would if we assumed the lHui fc Gnedinl (| 19971 ) 
rates, reducing the inferred clumping factor of ionised gas 
by the same ratio. The resulting uncertainty does not af- 
fect our results and is trivial to adjust for: Simply divide 
our inferred clumping factors by the ratio of the alternative 
recombination coefficient to our adopted coefficient at 10 4 
K. 

In summary, Figure indicates that the clumping fac- 
tor's value depends rather strongly on its definition. Taking 
self-shielding into account suppresses it from ~ 30 to « 6 at 
2 = 5 because the most overdense regions are self-shielded 
and do not contribute to the IGM recombination rate. Tak- 
ing the IGM's temperature into account further suppresses 
clumping to 2-4 because photoionised regions are generally 
warmer than 10 4 K. 



4 THE ROLE OF FEEDBACK 

Galactic outflows and photoheating both affect gas clump- 
ing. Heating suppresses gas clumping by smoothing den- 
sity flu ctuations on sca les smaller than the local Jeans 
length (|Efstathioul Il992h . Outflows increase IGM clump- 
ing partly by boosting the amount of dense gas that re- 
sides within halos but outside of galaxies and partly (for 
a given f eac {z)) by delaying reionisation and heating. The 
latter effect is a straightforward consequence of suppress- 
ing star formation. The former is expected based on our 
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Figure 4. The impact of photoheating and outflows on the over- 
all IGM baryonic clumping factor Cb. AH simulations subtend a 
6h _1 Mpc volume discretized with 2 X 256 3 particles (r6n256); see 
Tablc[T] Outflows have negligible impact on Cb in the absence of 
an EUVB, but when acting in concert with an EUVB their im- 
pact is stronger than uncertainties owing to limitations in our 
particular treatment of the EUVB. 

previous finding that outflows boost the nu mber of ionis- 
ing p hotons required to achieve reionisation (|Finlator et al.l 
l201lf ) although it may in reality be weaker if ejected gas 
self-shields. In this section, we quantify how these processes 
modulate gas clumping. 

4.1 Averaging over the IGM 

We begin by discussing the separate impacts of outflows and 
photoheating on the IGM clumping factor. In order to facil- 
itate this discussion, we revert to the definition Cb, which 
uses all gas that is too diffuse to form stars (Section[2]) with- 
out reference to its ionisation state or temperature. While 
this definition includes gas that is dense enough to self-shield 
and is therefore an imperfect estimate of the IGM recom- 
bination rate, it allows us to compare simulations with and 
without an EUVB. 

We use solid black and dotted blue curves in Figure|4]to 
show how Cb evolves in simulations without and with out- 
flows in the absence of an EUVB. They are nearly coincident, 
indicating that outflows affect the gas density distribution 
only weakly if they do not couple to an EUVB. Adopting 
a spatially-homogeneous lHaardt fc Madaul (|200ll . hereafter 
HM01) EUVB dramatically suppresses the clu mping factor 
follow ing 2 = 9 (long-dashed green; see also iPawlik et al.l 
120091) . 

Our self-consistent simulation omitting outflows 
(r6n256nWwRT16, short-dashed red) predicts a reionisa- 
tion history in which the neutral hydrogen fraction drops 
to 50% at 2 = 9.2, quite similar to the redshift at which the 
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HMOl EUVB reionises the Universe (z — 9). Prior to z = 5, 
its Cb exceeds the HMOl case owing to its more extended 
reionisation history. By z = 5, its Cb converges to the 
HMOl case although it shows several spikes that are likely 
associated with photoevaporation of star-forming gas from 
low-mass (M^ < 10 Mq) haloes in regions that are just 
being reionised. This suggests that Cb is weakly sensitive 
to the reionisation history at early times but converges to 
expectations from calculations using a spatially-uniform 
EUVB fairly rapidly (probably on the sound-crossing time 
of the filaments). 

Turning to our fiducial simulation (magenta dot- 
dashed), which does include outflows, we find that its Cb 
is higher than the no- wind case at all redshifts. This model 
has significantly less star formation than the no- wind model, 
but it also assumes a higher ionising escape fraction in order 
to compensate (Section 12. 2[) . Combining these factors leads 
to a reionisation history in which the neutral fraction drops 
below 50% at z = 8.9, only slightly later than the other 
models. Despite their similar reionisation histories, the wind 
model's Cb is roughly 30% higher at z = 5. Hence although 
the reionisation histories of the wind and no-wind models 
are not exactly the same, Figure [4] confirms that outflows 
enhance Cb by boosting the amount of gas that is overdense 
but not star-forming. 

We draw several conclusions from Figure [4] First, out- 
flows do not affect Cb unless they are heated by an EUVB. 
Intuitively, outflows consist largely of material entrained 
from their host galaxies' interstellar media, hence they re- 
main relatively cold and dense unless they are heated by an 
EUVB. Second, an EUVB heats and evaporates gas out of 
shallow potential wells such as filaments and minihalos, sup- 
pressing Cb- Finally, a better understanding of the nature 
of galactic outflows is at least as important to our under- 
standing of the IGM density structure as the nature of the 
EUVB. In Section f4. 21 we will ask whether this uncertainty 
impacts the recombination rate of diffuse gas. 
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Figure 5. (top) The gas density PDF in units of mass-weighted 
probability per unit log 10 (A) for six different simulations as indi- 
cated by the legend in the middle panel at z = 6; all simulations 
subtend volumes with side length 6/i -1 Mpc and resolve halos of 
total mass 1.4 X 10 s Mq. The vertical dashed line indicates the 
threshold density for star and wind formation (0.13cm -3 ). Pho- 
toionisation heating pressurizes gas with densities A > 10, and 
outflows transfer gas from galaxies into the IGM. (middle) The 
cumulative baryonic clumping factor as a function of maximum 
overde nsity compared to the fiducial simulation of iPawlik et al.l 
<|2009h (black dotted). An EUVB suppresses the clumping of all 
gas, with minor differences owing to the details of outflows and 
reionisation. (bottom) The neutral hydrogen fraction as a func- 
tion of overdensity. Gas above a threshold overdensity (which is 
related to Aj ls ) is neutral owing to our self-shielding model. 



4.2 The recombination rate as a function of 
density 

Figure U suggests that outflows and photoheating modulate 
the IGM's density structure, but their impact on the IGM re- 
combination rate remains unclear because Cb averages over 
gas that is in reality self-shielded. In order to understand 
how feedback processes modulate the abundance of gas at 
different densities, we show the mass-weighted probability 
density function of gas (PDF) at z = 6 as a function of over- 
density A = p/(p) in the top panel of Figure [5] Comparing 
the blue and black solid curves to the lower curves shows 
how photoheating smooths ove r filaments and rem oves gas 
from low-mass haloes (see also lPawlik et ai1l2009T ). At den- 
sities higher than the adopted threshold for star formation 
(vertical dashed line), the PDF is depleted by star formation 
and outflows. 

Above A ~ 240, both reionisation simulations (short- 
dashed red and dot-dashed magenta) show an increas- 
ing gas abundance because the EUVB is attenuated in 
these regions owing to self-shielding. Outflows modestly 
increase the gas abundance throughout the region 10 < 
A < 1000. The long-dashed green curve shows the PDF 
from the r6n256nWHM01 simulation, which omits outflows 



and assumes a spatially-homogeneous EUVB (with no self- 
shielding). Its gas abundance is suppressed with respect to 
the self-shielding calculations in regions where A > 1000 
and enhanced near A ~ 200 because the EUVB evaporates 
dense gas that should in reality be shelf-shielded. 

Differences in the gas density PDFs translate into dif- 
ferent amounts of gas clumping. In the middle panel of Fig- 
ure [S] we show the cumulative clumping factor as a func- 
tio n of the threshold ov erdensity C(< A t hr), defined follow- 
ing [Pa^hkeFal] (|2009h : 



C(< A thr ) = 



J Athr dAA 2 P(A) 



(0) 



dAV(A) 

Here, V(A) is the volume-weighted probability distribu- 
tion function of gas density. Note that each curve's value 
at the star formation threshold equals the corresponding 
simulation's Cb at z = 6 in Figure [4] The solid blue 
and black curves confirm that outflows have little impact 
on the clumping factor without an EUVB. Adopting a 
spatially- homogeneous EUVB (green long-dashed) repro- 
duces the IPawlik et all l|2009l 'l result (black dotted) as ex- 
pected, with slight differences likely owing to their addi- 
tion of extra heating around the reionisation redshift. Our 
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radiation-hydrodynamic simulations again predict very sim- 
ilar clumping factors. Outflows boost clumping as expected 
from the top panel, but the change is slight. In particular, 
we confirm the finding bv lPawlik et alj (|2009l ) that outflows 
change C(< A t hr) by < 10% when averaged over gas with 
A < 100. This is somewhat nontrivial given that our outflow 
prescription ejects roughly twice as much gas per unit stellar 
mass formed as theirs at this redshift. Overall, the impact 
of changing the EUVB and the outflow model on C(< Athr) 
is quite modest. The fact that the gas density PDF at low 
densities (A < 200) is relatively robust to these factors indi- 
cates that the recombination rate in the diffuse IGM — which 
is the important quantity for understanding the progress of 
reionisation — is only weakly sensitive to the uncertain de- 
tails of star formation and feedback. At higher densities, the 
impact of feedback and self-shielding is much stronger (top 
panel), but recombinations in these regions do not count to- 
ward the IGM recombination rate as their ionisation state 
is dominated by local sources. 

The ability of outflows to boost the abundance of neu- 
tral gas at moderate overdensities (A > 240) rather than 
only at the outskirts of galaxies (A > 1000) suggests that 
they could function as an additional population of Lyman 
limit systems residing near or within halos. If so, then they 
would have to be self-shielded against the EUVB. We show 
in the bottom panel the neutral hydrogen fraction as a func- 
tion of A for our fiducial simulation. The neutral fraction 
jumps from to 1 at A ~ 240, indicating that outflows 
could indeed be partially self-shielded. This possibility is 
consistent with direct observations of cold gas in gala ctic 
outflows (for example. lMartinll2005l; I Weiner et alj|200gh a s 
well as recent theoretical results ( van de Voort et al.l 2012T ). 

In summary, an EUVB significantly suppresses IGM 
clumping by pressurizing moderately overdense regions 
while outflows modestly boost the recombination rate by 
returning gas from galaxies into the IGM. Meanwhile, self- 
shielding hides the most overdense gas from the EUVB. This 
boosts the mass fraction of neutral gas while suppressing the 
clumping factor in ionised regions. Outflows are at least par- 
tially self-shielded, suggesting that they could contribute to 
the IGM opacity. 



5 APPLICATION I: PHOTON-COUNTING 

Analytical models of reionisation u se clumping factor s 
to account for IGM inhomogeneities (|Madau et al.l 1 1999T ) , 
prompting the need for theoretical insight into the nature of 
clumping. We have discussed the various effects that modu- 
late the value of the clumping factor that we derive from our 
simulations. In order to close the loop, we now ask whether 
these clumping factors do indeed constitute a good model 
for the IGM recombination rate during reionisation. We do 
this by using an analytical model to explore how well several 
different clumping factors fare in reproducing the reionisa- 
tion history of our r6n256wWwRT16 simulation. We will 
show that clumping factors perform reasonably well, and 
that they fare better if they incorporate information from 

the simulated temperature and ionisati on fields. 

We use as our model Equation 20 of lMadau et al.l (|l999l ) 
(this can also be derived by taking the spatial average of 
Equation [4] and normalizing by the mean hydrogen number 



density) : 



tlx 



HI,V 



dt 



+ recombinations. 



(10) 



Here, ihi,v = (ihi/wh) is the volume-averaged neutral hy- 
drogen fraction while h 1 is the ionising luminosity per hy- 
drogen atom into the IGM. We will use the value for n 7 that 
is predicted directly by our simulation. The form of the re- 
combination term varies depending on the clumping factor 
definition. 

There are three caveats to this widely-used formalism. 
The first is that the second term represents a mass-average 
while the third term is generally computed in a volume- 
averaged way (see, for example, Equation [SJ . To see this, 
recall that the second term is meant t o model the rate o f 
growth of the ionised volume fraction |Madau et alj|l999h . 
The rate at which an ionised volume Vi grows depends on 
the ratio of the ionising luminosity iV 7 to the gas density 
at the position of the ionisation fr o nt wg , Vi = N^/uh- 
In practice, however, iMadau et alj (|l999T ) substituted the 
volume-averaged gas density (tih) into the denominator. 
This approximation is only accurate if the density is ho- 
mogeneous or if the sources are widely separated (as in the 
case of helium reionisation; iMcQuinn et alj 120091 ). hence it 
is equivalent to assuming that the neutral mass and vol- 
ume fractions are equal, ihi.m = ihi,v- It is possible to 
use knowledge of the distribution of gas densities and the 
topol ogy of reionisation (for example, iMiralda-Escude et al.l 
120001 ) to compute the Vi from the ionising emissivity (or, 
alternatively, to compute the recombination rate from the 
ionised mass fraction), but these steps are omitted in an- 
alytical calculations. Owing to this approximation, we will 
refer only to the neutral fraction xhi rather than to Khi,m or 
ihi.v within the context of our analytical model. The sec- 
ond caveat is that Equation [TO] omits collisional ionisations 
because they are nontrivial to model analytically. Finally, 
the third term requires knowledge of the IGM's density and 
temperature distributions that is generally encapsulated via 
clumping factors. All of these simplifications are relaxed in 
our numerical calculations. In this section, we will use our 
simulations to ask whether they are indeed valid. 

We compare the resulting reionisation histories with our 
numerical simulation in Figure |S] As a baseline, the short- 
dashed blue curve illustrates the reionisation history if there 
are no recombinations (C = 0). Comparing it with the dot- 
dashed green curve indicates that recombinations delay the 
completion of reionisation from z « 10 to z ~ 7.5. The de- 
lay is expected given that this simulation requires 3.5 ionis- 
ing photons per hydrogen atom to reach a volume-averaged 
neutral fraction of xhi = 0.01. The other colored curves rep- 
resent three different ways of treating recombinations. The 
solid black curve uses the clumping factor of all IGM gas 
and computes the recombination term in Equation [10] as: 
Cbfci(10 4 K)a:Hii(nH), where cehii = 1 — ihi represents the 
ionised fraction. Down to z = 10, this crude treatment is 
already a significant improvement over the C = curve. 
Below z = 10, recombinations begin to win over the de- 
clining ionising emissivity, which is in turn driven by the 
assumed f csc (z); the result is a multimodal reionisation his- 
tory. The fact that the simulation does not, in fact, yield a 
multimodal reionisation history emphasizes the importance 
of allowing overdense regions to self-shield against the ion- 
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Figure 6. The dependence of the reionisation history on the def- 
inition of the clumping factor, as compared to our fiducial radi- 
ation hydrodynamic simulations. In photon-starved reionisation 
scenarios, a rapidly-growing clumping factor can lead to a doublc- 
reionisation history although this does not happen with our pre- 
ferred definition of the clumping factor. 



ising background. The dotted red and long-dashed magenta 
curves correspond to the definitions Chii and C HII10 4 K 
from Section [4.11 In these models, the recombination term 
is Cfei(10 4 K)a-Hn{wH). They are each a clear improvement 
over using Cb. Recall that the dotted red curve accounts for 
the order in which regions of different overdensities reionisc 
while the long-dashed magenta curve additionally accounts 
for the impact of temperature fluctuations. Including this 
information accelerates reionisation by suppressing recom- 
binations (Figure [3j). 

In detail, reionisation is more extended in the numeri- 
cal model even when the clumping factor accounts for both 
the temperature and the ionisation fields (long-dashed ma- 
genta). We attribute this to the fact that the analytical 
model effectively assumes an infinite speed of light c. We 
have previously shown that arti ficially boosting c ca uses 
reionisation to complete earlier l|Finlator et al.l l2009bh al- 
though the effect was slightly weaker (Az — 0.1) than what 
is seen here (Az m 0.9). This delay is only a few times larger 
than the light travel time across the photon me an free path 
at this epoch (~ 10 Mpc; [Furlanetto fc Oh|[^005l ) ■ consistent 
with the possibility that light-travel time effects delay the 
completion of reionisation. 

There are two additional effect that extend the final 
stages of reionisation in numerical models. The first results 
from the fact that the timescale for ionising an absorber 
varies inversely with its geometrical cross section. For ex- 
ample, an EUVB with amplitude T at the frequency corre- 
sponding to an absorption cross-section a ionises a spherical 
absorber of uniform density n and radius r on a timescale 
|ncrr/r compared to 1/r in the optically-thin case. Approx- 



imating Lyman limit systems as uniform spheres in hydro- 
static equilibrium at a temperature of 10 4 K, we find that 
systems with overdensities of 10 and 100 are ionised by a 
background T — 2 x 10 _13 s _1 on timescales of ~ 8 and 24 
Myr at z — 6, compared to ~ 0.2 Myr in the optically thin 
case. This effect is omitted in analytical models but present 
in numerical simulations, although it probably delays the 
completion of reionisation by only Az « 0.1. 

The second delay results from the fact that photons 
with high energy and small absorption cross section take 
longer to be absorbed than photons near the Lyman limit 
once the ionised volume fraction is large. This effect is miss- 
ing from our simulations owing to our monochromatic radi- 
ation transport solver. However, it is not large: for z 6, 
energies less than 54.4 eV, and neutral hydrogen fractions 
greater than 0.01, the delay is less than 20 Myr. 

We may quantify the error in the reconstructed reioni- 
sation history using the optical depth to Thomson scattering 
T es - We compute this for the simulations and the analytical 
models assuming that helium is singly- ionised with the same 
ionisation fraction as hydro gen down to z = 3 (|Ciardi et al.l 
|2012| ; iFriedrich et al.|[2012h . and doubly-ionised thereafter. 
In our simulation, the predicted value is 0.071. Ignoring re- 
combinations yields r cs = 0.095, which confirms that galax- 
ies could have dominated reionisation modulo uncertainties 
regarding the star formation efficiency of low-mass halos and 
the IGM recombination rate. Adopting realistic clumping 
factors measured directly from our simulation yields values 
for T cs that are between 0.062 (solid black) and 0.082 (long- 
dashed magenta). The discrepancy between the long-dashed 
magenta and dot-dashed green curves can be regarded as an 
estimate of the uncertainty in analytical reionisation calcula- 
tions owing to complications such as light-travel time effects, 
shadowing, and source clustering; it is w 10%. In short, an- 
alytical models of reionisation can reproduce the behavior 
of more realistic models reasonably well, with improvement 
possible if the adopted recombination rate accounts for the 
inhomogeneous temperature and ionisation fields. Our an- 
alytical modeling also confirms that collisional ionisations 
and the other complications that we have mentioned are 
subdominant. 

We note that iGnedinl |2000r i previously performed 
an extensive comparison between his numerical radiation 
hydrodynamic simulations and analytical models, finding 
broad support for the basic assumptions underlying ana- 
lytical models. Our results agree with his. The principal dif- 
ference is that, in computing the volume-averaged recombi- 
nation rate, their analysis did not distinguish between gas 
within galaxies and the IGM. This led to a much higher 
value for the clumping factor (100-200) than is appropri- 
ate for densities that are conventionally associated with the 
IGM (for example, A < 1000). Our simulations sidestep this 
complexity by using Au 8 as a physically-motivated model 
for the boundary between the IGM and the condensed gas 
whose absorptions are more appropriately attributed to / eS c- 



6 APPLICATION II: THE PHOTON COST OF 
REIONISING FILAMENTS 

While much has been learned abo ut the photo n cost 
of removing gas from minihalos (|Haiman et al.l l200ll ; 
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Figure 7. (bottom) The dependence of the overlap redshift on 
heating redshift: Delaying photohcating from z = 19.5 to z = 7.5 
delays overlap by Az = 3. (top) The reionisation photon budget 
per hydrogen atom as a function of heating redshift: Delaying 
photohcating can more than double the photon cost of reionisa- 
tion. The crosses correspond to a homogeneous IGM (C = 1). 



Barkana fc Loebll2002l ; IShapiro et aT]|2004l ; llliev et al.|[2005l : 

Ciardi et all 12000 ) . the ability of moderately-overdense gas 



to delay reionisation has received less attention. This is po- 
tentially important for two reasons. First, analytical mod- 
els may underestimate the reionisation photon budget by 
deriving the clumping factor from hydrodynamic simula- 
tions in which the IGM is effectively pre-heated. For exam- 
ple, it is possible to compute the clumping factor at high 
redshift by applying Equation [9] to t he baryonic density 



redshilt by applying liqu ation [9| to t he baryonic density 
PDF of iMiralda-Escude et all ((2000I ) (|Wvithe et all 12008 ; 
iPritchard et al.ll2010l ). However, that PDF derives from sim- 
ulations in which the diffuse IGM is already smoothed by 
an EUVB. In reality, there is a photon cost associated with 
smoothing the IGM because the initially cold IGM has a 
higher clumping factor (Figure [5]). Ignoring this cost may 
cause models to underestimate the photon budget of reion- 
isation. Alternatively, it is possible that the IGM is heated 
before reionisation is well under way. For example, an early 
X-ray background co uld suppress gas clum ping prior to the 
onset of reionisation (|Oh fc Haiman|[200l 'l. In this case, fil- 
aments would absorb fewer ionising photons than in models 
(such as ours) in which the primary heating source is UV 
photons. In this section, we estimate the magnitude of this 
uncertainty by asking how many ionising photons are re- 
quired to reionise filaments, or regions with overdensities of 
A =1-50. 

Our approach involves decoupling reionisation from 
heating and determining how reionisation changes if the 
universe is heated at different redshifts. If the universe is 
heated at high redshift (such as z — 19.5), then the photon 
cost of reionisation is smaller because the universe remains 



smooth after being heated. By contrast, if the universe is 
heated at a lower redshift (such as z = 7.5), then reionisa- 
tion is more expensive because more of its gas has condensed 
into filaments by the time that reionisation is under way. 
By comparing these histories, we gain some insight (albeit 
not completely self-consistently) into the photon price that 
the universe pays for waiting to ionise its filaments; this is 
probably comparable to the actual photon cost of reionising 
moderately-overdense regions. 

iPawlik et al.l |2009i ) provide fitting functions for the 
baryonic clumping factor in simulations where an instan- 
taneously imposed spatially-homogeneous EUVB is used to 
heat the IGM at a variety of redshifts. We have combined 
their fitting functions (choosing their C50) with the pre- 
dicted ionising emissivity from our r6n256wWwRT16 simu- 
lation within our analytical model and evaluated how reion- 
isation depends on the smoothing redshift. We show in 
the bottom panel of Figure [7] how the reionisation red- 
shift (where ihi 0.01) varies. Delaying smoothing from 
z — 19.5 to z = 7.5 delays reionisation by Az ~ 4.5. This rel- 
atively strong dependence owes to the extended reionisation 
history that is enforced by our choice of f esc (z). We have ver- 
ified that extracting our emissivity history from simulations 
that assume a constant / csc = 0.5 yields a smaller delay of 
Az = 1.1 (not shown) because the growth rate of collapsed 
matter (such as stars) is much faster than the growth rate 
of moderately overdense structures (such as filaments). 

The top panel maps the reionisation redshift into the 
number of ionising photons consumed per hydrogen atom. 
If the IGM is heated at z = 19.5, reionisation consumes 1.8 
photons per hydrogen. Observations suggest that reionisa- 
tion was well under way by z =9 -10 (|Pritchard et alj|20ld : 
iKuhlen fc Faucher-Giguerej|2012r ). Adopting this as the pho- 
toheating redshift yields a total photon cost of 2.5-2.9 pho- 
tons per hydrogen. Subtracting, we find that the universe 
pays a price of 0.7-1.1 photons per hydrogen for waiting 
to heat its filaments. Put differently, analytical models that 
assume that the universe is smoothed at a high redshift un- 
derestimate the total photon cost of reionisation by 0.7-1.1 
photons per hydrogen atom. 

These calculations may underestimate the total photon 
cost of reionisation because they do not account for gas that 
is condensed into minihalos, or halos with virial tempera- 
tures below 10 4 K (« 10 s M© for redshifts of 6-10). Prior 
to reionisation, much of the gas may have been condensed 
into halos with virial masses betw een 1 6 and 10 8 M Q . The 
clumping factors of IPawlik et al.l (|2009l ) were derived from 
simulations that resolved halos at the upper end of this mass 
range with 100 particles, hence lower-mass systems are ef- 
fectively missing. T hese halos can consu me up to 5 photons 
per hydrogen atom l|Shapiro et al.ll2004h . with the implica- 
tion that the amount of gas clumping could have been much 
higher if the IGM temperature did not exceed 1000 K prior 
to reionisation. 

On the other hand, an early X-ray background could 
deposit an early entr opy floor due to the large mean 
free path of X-rays ifOhl |200ll: IVenkatesan et alj 1200 ll : 
iRicotti fc Ostrikerll2004l : iMadau et alj |2004 ) ■ Such a back- 
ground is widely predicted in many models of reionisation 
and could arise from X-ray binaries, inverse Compton scat- 
tering in supernova remnants, or early min i-quasars. This 
process significantly reduces gas clumping (|Oh fc Haimanl 
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l2003h . and if it occurs well before the reionisation epoch 
then it could decrease the reionisation photon budget. 

In short, while filaments may have absorbed ~ 1 ad- 
ditional photon per hydrogen, the total number of photons 
required to complete reionisation remains uncertain owing 
to the unknown abundance of small absorbers. Future work 
incorporating X-ray heating and a wider dynamic range will 
be required in order to explore these processes. 



7 DISCUSSION 

Could galaxies have produced enough ionising photons to 
reionise the Universe? Could they have maintained the 
IGM's ionisation state at z = 6? Two major uncertainties 
hamper efforts to connect observations of galaxies with mea- 
surements of the neutral hydrogen fraction in order to ad- 
dress these questions. The first is the fraction of ionising pho- 
tons that escape galaxies' ISMs, / eac , and the second is the 
recombination rate of the IGM, which is often parametrized 
by a clumping factor C. The observed population of galaxies 
can maintain a n ionised IGM if the rat io C// eS c lies within 
the range 1-10 (jFinkelstein et alj|20l3 . Figure 15). 

We previously showed that our simulations reproduce 
the observed UV luminosity function of galaxies and self- 
consistentl y complete reionisatio n by z = 6 when we assume 
/ csc = 0.5 (|Finlator et alj|201lh . While encouraging, those 
simulations suffered from two drawbacks. First, they over- 
produced the observed EUVB amplitude at z ^ 6 and under- 
produced r es . Second, the radiation transport solver's lim- 
ited spatial resolution did not shield moderately-overdense 
regions from the EUVB, with the result that the predicted 
clumping factor of ionised gas was artificially high (w 12 
at z = 6). Consequently, completing reionisation (xhi,v < 
0.01) required up to 5 ionising photons per hydrogen atom. 

In this work, we have improved on those calculations by 
adopting a subgrid self-shielding treatment and a redshift- 
dependent / csc . Self-shielding limits the clumping factor and 
the reionisation photon budget by isolating overdense gas 
from the EUVB. Meanwhile, an evolving / cac (1) reconciles 
the model with the observed r cs by increasing the electron 
fraction at early times; (2) heats the IGM at an earlier red- 
shift, which suppresses gas clumping at later times; and (3) 
raises the Gunn-Peterson optical depth at z = 5 into im- 
proved agreement with observations. Together, we find that 
these modifications suppress C to values in the range 2.7- 
3.3 at z — 6 (the uncertainty comes from possibly overes- 
timating the overdensity at which gas fully self-shields; see 
below). For our adopted f csc (z), this leads to C// osc ~25- 
30 at z = 6, slightly too high for currently-observed galaxy 
populations to maintain a completely ionised IGM. Our sim- 
ulations overcome this barrier by resolving galaxies down to 
absolute magnitudes of Mtjv — — 15, which is 2-3 magni- 
tudes fainter than current limits for observing galaxies di- 
rectly but consistent w ith requirements i nferred from other 
star formation tracers jTrenti et al.ll2012j ). 

Despite these improvements, there are a number of ways 
in which we can build on this work in the future. First, we 
use the global EUVB to compute An s . Near sources, this un- 
derestimates the EUVB, leading to an overdensity threshold 
that is too low. In voids, it overestim ates the EUVB, lead - 
ing to too little self-shielding (see also lCrociani et akll201l[ ). 



It would be preferable to estimate the EUVB by averaging 
over regions that are a few times larger than the expected 
length scale of self-shielding systems in order to incorporate 
the EUVB's small-scale fluctuations into the ionisation field. 
This improvement would result in more-ionised overdensities 
and more-neutral voids (or in other words a more strongly 
inside-out reionisation topology), but the overall impact on 
C is difficult to predict. 

Second, our results would benefit from an improved 
understanding of the IGM's temperature T. The IGM re- 
combination rate is sensitive to T through the recombina- 
tion coefficient and through the minimum overdensity A 
at which gas is neutral, which in turn is proportional to 
the self-shielding threshold An s . For temperatures between 

10 3.5_ 1 q5.5 ] hi x T -0.77^ Meanwhile] C(< A thr ) OC A°' 25 

(Figure [5| and An a oc T ' 3 (Equation [2]), leading to a weak 
overall dependence of C(< A t h r ) oc T ' 075 . Combining these, 
we find that the overall IGM recombination rate varies as 

ji-0.7 

Uncertainty in T stems from the unknown latent heat 
of photoionisation and from resolution limitations. We have 
effectively assumed that reionisation heats gas to ~ 15, 000 
K. While this assumption yields a mean IGM temperature 
that is consistent with observations (Section [2j, the obser- 
vations are themselves uncertain owing to the necessity to 
correct the meas ured IGM temperat ure for the influence of 
nearby quasars |Bolton et al.l I2OI2T ) . Assuming a (harder, 
softer) ionising background would lead to a (hotter, colder) 
post-reionisation IGM. Additionally, the limited spatial res- 
olution of our radiation solver may cause our simulations 
to underestimate the post-reionisation gas temperature be- 
cause ionisation fronts are artificially broadened, increas- 
ing the amount of time during which partially-ionised gas 
can cool through collisional excit ation of neutral hydro- 
gen (|Miralda-Escude fc Reesl [r994l'). The magnitude of this 
effect is likely small compared to the uncertain latent heat 
of photoionisation. Nevertheless, improved measurements of 
the IGM temperature at z > 5 would constrain our model 
further. 

Third, our simulations expose somewhat too much over- 
dense gas to the EUVB, boosting the IGM recombina- 
tion rate. For example, the threshold overdensity for self- 
shielding at 10 4 K and z = 6 should be An B = 25 (Equa- 
tion [2]). If the neutral column density of self-shielded sys- 
tems varies strongly with A, then systems with overdensity 
more than a few times greater than An a should be essen- 
ti ally neutral. For examp le, the high-resolution calculations 
of lMcQ uinn et ail |201lT ) indicate that the neutral fraction 
should increase rapidly above A ~ 100 (their Figure 3). By 
contrast, the gas in our simulations remains ionised until 
A ~ 200 (Figure [5]), which may be slightly too high. If so, 
then it owes largely to the fact that our simulations assume 
that all photons have an energy of 1.3 times the threshold 
energy for photoionisation /i^ll in order to photoheat the 
IGM. As we argued in the discussion of Figure [2l this as- 
sumption may artificially boost An s by a factor of 1.7. Invok- 
ing once again the scaling C oc Af ls 25 (Figure[S]), we find that, 
if the amplitude of the ionising background were unchanged 
but the radiation field were dominated by photons with en- 
ergy nearer to /i^ll, then the clumping factor of ionised 
gas would be lower by « 12%. Combining this uncertainty 
with the range of values inferred from our fiducial and high- 
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resolution simulations (r6wWwRT16 and r6wWwRT32), we 
conclude that the clumping factor at z = 6 lies between 
2.7-3.3. 

Finally, our limited cosmological volume forces us to 
adopt a rather dramatic evolution for / csc (z). We do this 
because our predicted r cs and r a are not converged with 
respect to simulation volume or mass resolu tion, as is well 
known for simulation volum es of this size (|Gnedin fc Fan! 
l200d : iBolton fc Becker! 120091 ). This is not to say that the 
need for a decreasing f csc (z) is purely a resolution limita- 
tion; empirical arguments strongly support such a trend and 
in fact our ado pted relation lies within the obser vationally- 
inferred range (jKuhlen fc Faucher- Giguerd[2012T ). However, 
with a larger volume, structure formation would begin 
sooner, allowing us to adopt a weaker redshift-dependence 
for /e S c (2) without compromising the agreement between the 
observed and predicted r ca and r a . This would change the 
predicted IGM temperature-density relationship as well as 
the overall gas density distribution. Both of these changes 
could modify the predicted IGM recombination rate, al- 
though the sign of the effect is difficult to predict. 



8 SUMMARY 

We have used a suite of cosmological radiation hydrody- 
namic simulations to explore the impact of outflows and 
an EUVB on the recombination rate of the reionisation- 
epoch IGM. Feedback processes modify the IGM's gas den- 
sity, ionisation, and temperature fields. We may illustrate 
their impact on the gas density distribution by considering 
the clumping factor C100 = A 2- P(A) at redshift z = 6 and 
overdensity A = 100; the overall clumping factor Cb (which 
neglects the ionisation and temperature fields) scales with 
this value. 

• In the absence of feedback, baryons follow dark mat- 
ter, a nd C100 = 10.7 a s previously found in N-bo dy simula- 
tions (|lliev et al.ll2007l ; iRaicevic fc TheundkOllT ). 

• Galactic outflows return star-forming gas into the IGM, 
boosting Cioo- The increase is less than 30% for any EUVB 
amplitude (including zero). 

• An ionising background pressurizes the diffuse IGM, 
suppressing Cioo to 2.0. 

In reality, the IGM's volume-averaged recombination rate 
depends on its density, ionisation, and temperature fields. 
Taking this information into account one piece at a time 
within the context of our fiducial simulation, we find (still 
at z = 6): 

• The clumping factor of all baryons below our adopted 
threshold density for star formation is Cb = 16.4. 

• Averaging only over the ionised gas yields a much lower 
value of Ch ii = 4.9 because overdense gas is self-shielded. 

• Accounting for the temperature field further suppresses 
the recombination rate to C HII 10 4 K = 3.3 because gas in the 
recent aftermath of reionisation is somewhat hotter than the 
canonical 10 4 K. 

• The temperature-corrected clumping factor of ionised 
regions increases from 0.79 at z = 15 to 3.3 at z = 6. Equa- 
tion [8] reproduces this evolution reasonably well throughout 
z = 15 -»• 5. 



• Our fiducial simulation may overestimate the minimum 
density of neutral gas. Correcting for this may suppress 
C HII 10 4 K by an additional factor of 12%. Hence our most re- 
alistic estimate of the ionisation- and temperature-corrected 
clumping factor at z — 6 is that it lies within the range 2.7- 
3.3. 

We have constructed an analytical reionisation model 
and tested how well different definitions of the clumping 
factor reproduce our numerical simulation's reionisation his- 
tory. We find that the clumping factor averaged over all 
IGM baryons Cb significantly overestimates the IGM recom- 
bination rate because it treats self-shielded gas as if it were 
optically thin. Accounting for the ionisation field through 
Ch 11 improves agreement with the numerical simulation, al- 
though the recombination rate is still overestimated because 
the post-reionisation IGM is in reality slightly hotter than 
the canonical 10 4 K. By contrast, using C HII10 4 K to ac- 
count additionally for the IGM temperature field yields a 
monotonic reionisation history whose r es overestimates the 
numerical result by only 10%. We speculate that this small 
remaining difference owes at least partly to light-travel time 
effects that are missing from the analytical calculation al- 
though other effects may also contribute. 

We have used our analytical model to estimate the pho- 
ton cost of reionising filaments, or regions with A =1-50. 
With our simulated emissivity history, if some process (other 
than reionisation) heats the IGM at z = 19.5 then the pho- 
ton cost of reionisation is y/H — 1.8 photons per hydrogen 
atom. Delaying the redshift at which the IGM is heated from 
z = 19.5 to z = 9.0 increases y/H to 2.9 because the clump- 
ing factor increases. The difference between these numbers 
indicates that the Universe pays a price of ~ 1 photon per 
hydrogen in order to reionise its filaments. 

The ionisation state of dense gas has observational im- 
plications that merit further study. Our simulations show 
that, in the presence of self-shielding, outflows can boost 
the amount of gas at densities corresponding to the cir- 
cumgalactic medium (10 < A < 1000). If a signifi- 
cant fraction of this gas remains neutral, then it could 
play two important roles in modulating the growth of 
the EUVB. Within halos, it could constitute a signifi- 
cant absorbing column, modifying the fraction of ionis- 
ing photons that escape the host halo. This effect could 
mimic a redshift-d ependent / ea c , which seems to be required 
by observations (llnoue et alj 120061; iFinlator et aL I l201ll ; 
iHaardt fc Madaul |2012| ; iKuhlen fc Faucher-Giguerd 120121 ). 
Additionally, gas that travels past the virial radius while re- 
taining a self-shielded component could b oost the abundance 
of Lym an limit systems, as suggested by Ivan de Voort et alj 
(|2012f ). 

A complementary motivation for improving our under- 
standing of self-shielded gas is that it can be observed di- 
rectly in the form of damped Lyman a systems and low- 
ionisation metal absorbers. The IGM's ionisation field is 
a major theoretical uncertainty hampering efforts to inter- 
prets ^obs^rvaH£ns_of_low-iojrisation metal ions. For exam- 
ple, (o^nenheimer^etji!] (|2009T l studied low-ionisation metal 
absorbers through the use of homogeneous ionising back- 
grounds. While this approach may be appropriate for ions 
such as SilV and CIV, it is less appropriate for low-ionisation 
species such as OI and Sill because it artificially ionises, 
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heats, and smooths gas that should in reality be self- 
shielded. Our current simulations allow dense gas to self- 
shield, hence they will yield improved predictions for the 
abundance of low-ionisation metal absorbers. This will re- 
duce the theoretical parameter space and increase the power 
of current observations. 
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Table 1. Our suite of simulations. The fiducial simulation is indicated in bold. 



name 


L a 


winds 


RT grid 


M h,min/ M O b 


r6n256nWnRT 


6 


no 




1.4 x 10 s 


r6n256nWwRT16 


6 


no 


16 3 


1.4 x 10 s 


r6n256wWnRT 


6 


yes 




1.4 x 10 8 


r6n256wWwRT16 


6 


yes 


16 3 


1.4 x 10 s 


r6n256wWwRT32 


6 


yes 


32 :i 


1.4 x 10 s 


r6nl28wWwRT16 


6 


yes 


16 3 


1.1 x 10 9 


r3nl28wWwRT8 


3 


yes 


8 3 


1.4 x 10 8 


r6n256nWHM01 


6 


no 




1.4 x 10* 



a in comoving h 1 Mpc 

b virial mass of a halo with 100 dark matter and SPH particles. 



APPENDIX A: APPROXIMATIONS TO THE EDDINGTON TENSOR 

Our method of solving the radiative transport equation involves using a time-independent ray-tracing calculation to update 
the Eddingto n tensors in a way tha t tracks the evolving emissivity and opacity fields without becoming too computationally 
burdensome (IFinlator et al.1 2009ah . The radiation transfer simulations in Table [1] use two approximations beyond those 



described in iFinlator et all (|2009aT ) 



First, we smooth each element of the Eddington tensor field with a 27-cell cubical tophat filter. This is necessary because 
the Eddington tensor can change rapidly at the positions of io nisation fronts, which c an in turn cause numerical errors when 
the radiation field is updated. We introduced this approach in IFinlator et"afl (|2009al ) but have not previously used it within 
the context of cosmological simulations. We stress that it does not materially degrade our spatial resolution because the 
smoothing occurs over the same length scale as the finite-difference stencil that we use to discretize the moments of the 
radiation transport equation. 

Second, we update each cell's Eddington tensor less frequently as it becomes optically thin. Here, the assumption is that 
a region's Eddington tensor evolves rapidly as an ionisation front sweeps over it and more slowly afterwards. We implement 
this idea as follows: The code only updates a cell's Eddington tensor if the cell's photon number density changes by more 
than a factor fj. Whereas previously we set fj = 0.05, we now allow fj to depend on the cell's optical depth r = xAirt, 
where \ an d Ax-rt are the local opacity and the cell size, respectively. In particular, fj — 0.05 if r ^ 1 and it grows linearly 
from 0.05 to 1 as r decreases from 1 to 0: 

_ / 0.05 t>1 (AU 

fj -\ 1.0 -0.957- r<l (A1) 

With this assumption, the Eddington tensor is updated frequently before the cell is reionised and less frequently once the 
radiation field has been established. 



APPENDIX B: RESOLUTION LIMITATIONS 

In this section, we assess the sensitivity of our predicted clumping factors to three different flavors of resolution limitation: 
The simulation volume, baryonic mass resolution, and the spatial resolution of our radiation transport solver. We evaluate 
convergence in Figure |BT1 by comparing the baryonic and observational temperature-corrected clumping factors from different 
simulations with different numerical parameters but the same physical assumptions. 



Bl Simulation Volume 

Our limited simulation volume coul d impact our results b y modifying the relative c ontribution of voids and overdensi- 
ties t o the overall clumping factor (|Bolton fc Beckerl 120091 : iRaicevic fc Theunsl l201ll ). By delaying the growth of struc- 
ture jBarkana fc Loeblliooil , small volumes also change the dependence of ionisation and temperature on overdensity at 
fixed redshift. Note that we have mitigated some of these problems by tuning f esc (z) so that our simulations match a variety 
of observations (although n ot all observations; S ection [2]), with the result that the simulations are representative even if they 
are not converged (see also lGnedin fc Fan]|2006t ). Nonetheless, it is of interest to test our sensitivity to our limited volume. 
The solid magenta and long-dashed green curves in Figure iBll correspond to simulations that subtend different cosmological 
volumes at the same physical resolution, and their predicted temperature-corrected clumping factors of ionised gas (top panel) 
are reasonably close. In detail, the simulation with the larger volume (solid magenta) predicts a slightly higher recombination 
rate at all times, with the gap growing to 20-30% by z = 5. This difference owes entirely to the fact that, in a larger volume, 
reionisation begins sooner, hence at any given redshift it has proceeded farther from the voids into the overdense regions. 
This is especially true in our simulations because our self-shielding treatment effectively enforces an outside-in reionisation 
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20 15 10 5 

redshift 

Figure Bl. The impact of various resolution limitations on the temperature-corrected clumping factor of ionised gas. The simulation 
volumes' side lengths L are in h~ Mpc, the SPH particle masses arc in 10 5 Mq, and the radiation transport grid sizes are in comoving 
h~ ^pc. Resolution limitations are generally weak, indicating that our fiducial simulation volume (solid magenta) is reasonably converged 
numerically. 



topology on small scales. In other words, for the same dependence of ionisation fraction and temperature on overdensity, the 
predicted recombination rates would be converged. 

Turning to the bottom panel, we see that the small- volume simulation has a higher volume- averaged recombination rate 
Cb- This may also be attributed to the delayed progress of reionisation: At any given redshift, photoheating has smoothed 
the density fluctuations up to a lower density threshold in the small-volume simulation because the ionising background is 
weaker. Consequently, the clumping factor averaged over all IGM baryons (including the self-shielded ones) is higher. 

This test does not account for fluctuations in the density field that are on larger scales than our larg est simulation volume , 
which is comparable to the typical size scale of ionised regions when the Universe is 50% ionised (|Furlanetto fc Oh 2005). 
Similarly, it does not account for fl uctuations in the EU VB on large (> 10/i _1 Mpc) scales, which can in turn source large-scale 
fluctuations in the ionisation field (|Crociani et aljfeoill ) . Future work at higher dynamic range will be required to assess these 
limitations. 



B2 Baryonic Mass Resolution 

The mass resolution of our SPH discretization dictates the amount of dense structure that can form, which in turn impacts the 
recombination rate by regulating the progress of reionisation. The solid magenta and short-dashed blue curves in the top panel 
of Figure IB II illustrate these effects. The result of reducing the mass resolution at constant volume is qualitatively similar to 
reducing the simulation volume at constant resolution because both delay reionisation. In particular, the curve corresponding 
to the higher-resolution simulation (solid magenta) climbs systematically above its lower-resolution counterpart below z = 13 
because reionisation has proceeded into regions of higher overdensity at a given redshift, boosting the mean recombination 
rate. Following reionisation (z < 7), the recombination rate in the lower-resolution simulation rejoins the prediction from 
the high-resolution simulation because the predicted gas density distributions and ionising backgrounds are similar. Turning 
to the clumping factor of all baryons Cb in the bottom panel, we see that the simulation with lower mass resolution has 
lower clumping prior to reionisation because the amount of matter that collapses into filaments and halos is lower at lower 
mass resolution. Following reionisation, Cb is higher at lower mass resolution because its weaker ionising background does not 
penetrate as far into overdense regions. 

It is important to note that varying the baryonic mass resolution changes the predicted clumping even if the spatial 
resolution of the radiation solver is unchanged. This illustrates the benefit of solving the cooling and ionisation equations on 
the SPH particles rather than on the RT grid. On the other hand, the change is not large despite the order of magnitude 
increase in resolution, indicating that the requirements for numerically resolving the clumping factor are not strict and that 
our simulations are therefore reasonably converged. 
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B3 Radiation Field Resolution 

A central benefit of modeling a spatially-resolved EUVB is the ability to resolve self-shielding within overdense regions. Once 
reionisation is well under way, the boundaries of shelf-shielded regions dominate the IGM opacity and the clumping factor. 
Increasing the radiation field's resolution treats self-shielding more accurately, preventing overdense regions from becoming 
photoionised and smoothed. While our subgrid self-shielding prescription prevents the most overdense regions from being 
photoionised, which limits the clumping factor, our radiation solver's grid size is sufficiently large that it could still ionise 
regions near the self-shielding threshold too efficiently; the result would overestimate gas clumping. Additionally, modeling the 
radiation field with a low spatial resolution artificially broadens ionisation fronts, extending the time during which partially- 
ionised gas cools through collisional excitation of neutral hydrogen and underestimating the post-reionisation temperature. 
Our fiducial simulation's spatial r esolution corresponds to an op tical depth of 2 at the Lyman limit and the mean density, 
hence this could affect our results |Miralda-Escude fc Reeslll994l ) In order to explore these possibilities, we compare the solid 
magenta and dotted red curves in Figurc lBTl which show how the clumping factor evolves when the radiation field is discretized 
using our fiducial and twice the fiducial resolutions, respectively. The predictions are nearly indistinguishable, indicating that 
the clumping factors are insensitive to the resolution of our radiation transport solver. This resolution convergence is a clear 
demonstration of the power of our subgrid self-shielding prescription. 

There are other numerical effects that we cannot consider directly using convergence tests. First, our use of a monochro- 
matic radiation solver likely results in underestimat ing the IGM temperature, particularly in voids, because it does not 
capture spectral hardenin g |Abel fc Haehneitl I1999T ) or heating owing to absorption of high-energy photons by helium 
atoms (|Ciardi et alj l2012h . This would slightly suppress the IGM recombination rate. On the other hand, the inability 
of ou r hydrodynamic sol v er to resolve minih a los means that o ur simulations ma y underestimate the IGM recombination 
rate (|Haiman et all [200H ; IShapiro et al.l |2003 ; llliev et all 120071 ; ICiardi et al.l 120061 ) ■ Resolving the question of whether nu- 
merical limitations cause us to over- or underestimate the clumping factor will require simulations with significantly higher 
dynamic range. 



